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ABSTRACT 

We report the results of 1-D hydrodynamical modelling of the evolution of gas in galaxy clus- 
ters. We have incorporated many of the effects missing from earlier 1-D treatments: improved 
modelling of the dark matter and galaxy distributions, cosmologically realistic evolution of the 
cluster potential, and the effects of a multiphase cooling flow. The model utilises a fairly standard 
1-D Lagrangian hydrodynamical code to calculate the evolution of the intracluster gas. This is 
coupled to a theoretical model for the growth of dark matter density perturbations. The main 
advantages of this treatment over 3-D codes are (1) improved spatial resolution within the cooling 
flow region, (2) much faster execution time, allowing a fuller exploration of parameter space, and 
(3) the inclusion of additional physics. 

In the present paper, we explore the development of infall models - in which gas relaxes into a 
deepening potential well - covering a wide range of cluster mass scales. We find that such simple 
models reproduce many of the global properties of observed clusters. Very strong cooling flows 
develop in these 1-D cluster models. In practice, disruption by major mergers probably reduces 
the cooling rate in most clusters. The models overpredict the gas fraction in low mass systems, 
indicating the need for additional physical processes, such as preheating or galaxy winds, which 
become important on small mass scales. 

1 INTRODUCTION 



Rich clusters of galaxies have long been known to be strong sources of diffuse X-ray emission, caused 
by thermal bremsstrahlung from a hot (T ~ lO^K) intracluster plasma (the intracluster medium, 
or ICM). This gas is the dominant baryonic component in such systems, contributing 20-30% of 
the gravitating material within an Abell radius. More recently, it has become clear that such a 
diffuse component is also present in poor clusters and galaxy groups. A good understanding of the 
processes involved in the formation and evolution of this diffuse gas is therefore of key importance 
for any complete theory of structure formation over a large range of scales. 

The earliest numerical models of the evolution of the intracluster gas (Gull & Northover 1975; 
Lea 1976; Takahara et al. 1976) were all one-dimensional simulations, concerned primarily with 
the evolution of uniform density gas 'clouds' relaxing into static potentials formed by the cluster 
galaxies and dark matter. In this simple situation, infall of gas leads to a rise in central gas density, 
and a shock propagates outwards from the centre of the cluster, heating the gas to temperatures 
similar to those observed today. After passage of the shock, the gas is approximately in hydrostatic 
equilibrium and further evolution is quasistatic. The subsequent cluster luminosity is found to be 
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nearly constant, although Lea (1976) found large fluctuations in Lx, probably due to numerical 
problems. The resulting gas distributions are generally not well fitted by any polytropic distribu- 
tion, although Gull & Northover (1975), who considered large radius infall into an initially empty 
potential, found that the cluster gas distributions were almost adiabatic. 

Perrenod (1978) considered more realistic infall models in which the cluster potential varies 
with time, using the results of an A^-body simulation of the evolution of the galaxies in a Coma- 
like cluster (White 1976). In these simulations, the cluster galaxies collapse rapidly by violent 
relaxation, followed by a slower contraction due to two-body encounters. Unlike models with a 
static potential, the cluster X-ray luminosity was found to rise by an order of magnitude from 
z = 1 to z = (for a flat Einstein-de Sitter cosmology), due primarily to the deepening of the 
potential well. However, this rate of evolution is considerably larger than that observed and may 
reflect ovcrmcrging in the galaxy distribution - White's models assumed that all the cluster mass 
was bound to the galaxies, overestimating the strength of two-body encounters. 

Hirayama (1978) extended the purely primordial simulations and considered models including 
mass and energy injection from a static galaxy distribution. These simulations showed that the 
injected material dominated the ICM towards the centre, primarily because the gas distribution 
in the cluster was more extended than that assumed for the galaxies. For high gas injection 
temperatures {T^nj ^ 3 x 10^ K) a cluster wind formed. 

More recently, much effort has been expended on the development of 3-D codes that are able 
to calculate the evolution of the gas and one or more collisionlcss components (e.g. dark matter, 
galaxies) from the primordial density fluctuation spectrum. This allows a fully self-consistent treat- 
ment in which the ICM evolves in a three-dimensional, time-varying potential. Since hierarchical 
clustering models will, for many interesting power spectra, lead to a high degree of substructure 
at all stages of the evolution, grid-based hydrodynamical techniques with their limited spatial 
resolution are unable to evolve the different levels of clustering hierarchy simultaneously. The gas 
flow is therefore generally calculated using particle-based techniques, for example smoothed parti- 
cle hydrodynamics (SPH; Gingold k. Monaghan 1977), coupled with a suitable AT-body algorithm 
(e.g., P^M; Efstathiou k Eastwood, 1981). 

Several authors have computed 3-D evolutionary models using such techniques - e.g., Evrard 
(1988, 1990), Thomas & Couchman (1992), Katz & White (1993), Metzler & Evrard (1994), 
Navarro, Prenk &; White (1996), hereafter NFW96. In most cases, cluster cores are found to 
develop at the intersection of filaments and sheets, and cluster growth proceeds through accretion 
of surrounding material as well as through mergers with smaller, collapsed systems. 

Although the general problem of non-linear evolution requires resort to numerical techniques, 
it is possible to study the evolution of structure in hierarchical models of the Universe by exploiting 
the approximate self-similarity of the evolution expected under certain conditions: namely those 
in which the background universe has closure density (fi = 1), so the Universe is scale free, and in 
which the amplitude of the initial density fluctuations follows a power-law, so these perturbations 
are also scale free. In this situation, the clustering is self-similar in the sense that statistical 
properties of the cluster population must evolve according to simple scaling laws (Kaiser 1986). 

In the idealised case of the spherically symmetric growth of an overdense region within an 
environment in which the overdensity varies initially as a power law function of radius, the evo- 
lution of individual clusters develops in a self-similar fashion (Fillmore k, Goldreich 1984 (FG84); 
Bertschinger 1985). In this case, the profiles of a growing cluster at different epochs are similar, 
and clusters of different mass at a given epoch arc simply scaled versions of one another. Since it 
has been shown by Hoffman k Shaham (1985) that the mean primordial overdensity distribution 
about a local peak does take a power law form, one might expect such self-similar growth to give 
a reasonable approximation to the mean growth of a cluster. 

In practice, the similarity will be disturbed by departures from spherical symmetry, and by 
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interaction between neighbouring density peaks. Nonetheless, simulations (e.g. Navarro, Frenk 
& White 1995, NFW96) indicate that when additional physics, such as radiative cooling and 
injection of winds from galaxies, is not included, cluster growth is approximately self-similar, with 
density profiles similar in form to those predicted by infall models of cluster formation (FG84; 
Bertschinger 1985). Also, clusters of different mass are very similar when appropriately scaled - 
though when a wide mass range is considered, lower mass systems are found to be rather more 
centrally concentrated, since they tend to form earlier. These results encourage the idea that 1- 
dimensional simulations should give a reasonable approximation to the mean evolution of clusters. 

Since previous 1-D treatments were too simplistic (none combined a realistic evolving potential 
with gas injection and cooling), and the limits of present day computers mean that 3-D simulations 
have poor resolution ( > 100 kpc), numerical simulations have thus far been unable to shed much 
light on a number of unanswered questions: e.g.. Do the scaling laws break down when cooling is 
included? Were cooling flows more massive in the past? How does the mass deposition rate vary 
with radius and system size? How does the gas mass fraction vary with radius? How does the 
cluster luminosity vary with time? 

The purpose of the present study is to address some of these questions, using a 1-D code, 
EVOL, that incorporates many of the effects missing from the older 1-D treatments: improved 
modelling of the dark matter and galaxy distributions, more realistic gas injection mechanisms, 
the effects of convection and multiphase cooling, as well as a simple investigation of the effects of 
mergers. 

The model utilises a fairly standard 1-D Lagrangian hydro dynamical code to calculate the 
evolution of the intracluster gas. This code is coupled to a theoretical model for the growth of 
dark matter density perturbations, which gives a cluster potential that deepens in a cosmologically 
realistic way with time. The main advantages of this treatment over 3-D codes are (1) improved 
spatial resolution within the cooling flow region, (2) fast execution time, allowing a fuller explo- 
ration of parameter space, and (3) the inclusion of additional physics. The main disadvantage is 
that we are unable to accurately model the effects of asymmetric processes, such as cluster merg- 
ers. This limitation, however, does not invalidate our results since, although the evolution of the 
dark matter in clusters is complex, with many clusters showing evidence of substructure, the dark 
matter distribution in "quiescent" clusters can be reasonably well approximated by a spherically 
symmetric model. For example, recent simulations by Tormen, Bouchet & White (1996) show 
that dark matter halos spend about two thirds of their evolution in a relaxed state, during which 
the density profiles of the halos are well approximated by the NFW96 analytical model. Having 
said this, periodic major mergers may well have a substantial impact on the structure within the 
central cooling flow region, possibly disrupting the flow. We investigate this effect below, using a 
simple model for the change in entropy of the gas within the cooling flow following a merger. 

Our aim is to determine the role of, and place constraints on, the heating and cooling processes 
that are important for the evolution of the ICM. The current paper deals with models in which the 
intracluster gas is unaffected by injection of energy and heavy elements from cluster galaxies, in 
order to test the extent to which the observed properties of clusters and groups can be explained by 
models with gas dynamics and radiative cooling, without recourse to additional physics. Models 
with gas injection from the galaxies are considered in a second paper (Paper II - Ponman and 
Knight, in preparation). 

The layout of the paper is as follows: the equations on which the code is based, and the 
numerical methods employed, are outlined in Section ^, and some test runs used to verify the 
accuracy of the code are described in Section ^. A large number of simulations have been run, 
in order to investigate the mass dependence of the evolution, and to check the effect of varying a 
number of key physical parameters. This grid of models is described in Section ^, and the evolution 
of the hot gas which results is discussed in broad terms in Section Some more detailed aspects 
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of the properties of these simulated clusters are then investigated: cooling flows in Section y, the 
distribution of hot gas in Section and the evolution of the X-ray luminosity in Section |8|. The 
simulations presented here assume anQ = 1 Universe, however in Section ^ we examine briefly the 
impact of a lower value of 0, on the X-ray evolution, and conclude that it is not large. Finally, in 
Section IG, we summarise our conclusions. 



2 NUMERICAL METHOD 

2.1 Basic Equations 

Assuming spherical geometry, the evolution of the ICM in galaxy clusters is described by the 
standard time-dependent equations of gas dynamics (Mathews & Baker 1971), with additional 
terms to remove gas from the flow when the cooling time is sufficiently short in order to simulate 
the features of an inhomogeneous cooling flow. Under such conditions, the equations of mass, 
momentum, and energy conservation that govern the evolution of the ICM are 

dp 1 d(r^pu) 



Dt p dr 
De p + q Dp rienn 



A, (3) 



Dt p2 JJf p 

where D/Dt is the Lagrangian derivative, u the radial gas velocity (positive for outflow), p the ICM 
mass density, Ue the electron number density, nn the hydrogen number density, and lo accounts for 
gas lost by cooling (see section 2.2 below). For a gas with ratio of specific heats 7 = 5/3, the specific 
energy e and the gas pressure p are equal to 3kT/{2fimp) and pkT/{fj,mp) respectively, where T 
is the gas temperature, k is the Boltzmann constant and nip is the proton mass. Throughout this 
paper we take // = 0.6 and /Xe = 1.15, appropriate for a fully-ionized plasma composed of 90% 
hydrogen and 10% helium by number. 

The term q in Equations § and |^ is the artificial viscosity (Richtmyer & Morton 1967), 

q = p(amml.^,o}] , (4) 



dr 

which is effective at smoothing the boundaries of shock waves. The parameter a determines the 
thickness of the shock front. Setting a = 2, the value that is adopted throughout this paper, 
ensures accurate modelling of the shock by spreading the jump in parameter values over several 
Lagrangian shells. 



2.2 Radiative Cooling and Mass Deposition 

The term on the right hand side of the energy equation accounts for the effect of radiative cool- 
ing. A is evaluated using bilinear interpolation from tabulated versions of the cooling function of 
Raymond, Cox, &: Smith (1976) for primordial and solar abundances, 

A{T,Z} = {1-Z)A{T,0} + ZA{T,1}, (5) 

where Z is the metallicity of the gas. Compton cooling is only important at high redshift (z > 10), 
and is ignored in the simulations presented here. At such early epochs additional energy input 
from early star formation is likely to be effective at preheating the cluster gas. The effects of this 
will be explored in Paper II. 



4 



Observed X-ray surface brightness profiles in cooling flows imply that gas cools and is deposited 
throughout the flow (Fabian, Nulsen &: Canizares 1991). The treatment of mass deposition used 
in EVOL is similar to that of David, Forman & Jones (1990). Gas is assumed to be thermally 
unstable to entropy perturbations when (2 — dlnA/dlnT) > (David, Forman & Jones 1990). 
When this condition is satisfied, and the integrated isobaric cooling time, tcooi is less than a^dt, 
gas is removed from the ICM at the rate Cp/trad, where p is the gas density, tj-ad is the local 
instantaneous cooling time, and Omd and ^ are model parameters. Theoretical considerations and 
observations of the cooling flows in rich clusters suggest ~ 1 (Sarazin 1990). This accounts for 
the last term in Equation ||, 



where M^arkii"), M[CM{r) and M^epif) are the total dark matter, ICM and 'cooled' gas mass 
within radius r respectively. The cooled gas mass, M^^p^ is primarily used for book-keeping, 
keeping track of material that is lost from the ICM due to mass deposition, although in the 
simulations it was found to be dynamically important close to the cluster centre at later times. 
The dark matter distribution otherwise dominates the dynamical evolution of the ICM. The galaxy 
distribution is only important in models with significant feedback from galaxy winds and ram- 
pressure stripping and is not considered here. (The galaxies are effectively incorporated into the 
dark matter component.) 

The evolution of the dark matter in clusters is complex, with many clusters showing evidence 
for recent mergers. However, the dark matter distribution in quiescent clusters is reasonably 
well approximated by a spherically symmetric model. NFW96 found that the radial dark matter 
density profiles of their simulated clusters were approximately self-similar, with profiles at large 
radii similar in form to that predicted by analytic infall models (FG84; Bertschinger 1985). In these 
models, the physical mass distribution of a cluster growing by the 1-D accretion of surrounding 
material at time t may be calculated by scaling a canonical function to the characteristic size 
{Rvirit)) of the system. We assume that such models give a good approximation to the shape and 
evolution of the dark matter distribution of real clusters, and adopt the canonical profile given by 
the similarity solution of FG84 corresponding to n = —1 (see Hoffman & Shaham 1985), which 
matches the slope of the CDM power spectrum on cluster scales (Kaiser 1986). This model tends 
towards the critical background density pc at large radii. Cusps in the density arising from the 
turnaround of particle shells are present in the similarity solution. In practice, such features would 
be washed out by departures from perfect spherical symmetry - we therefore eliminate them by 
fitting a power law profile to the formal similarity solution. 

Following the n = —\ similarity solution, the full dark matter profile at z = is determined by 
a single parameter, which may be taken to be the total mass inside the virial radius R^ir (defined 
throughout this work as the radius within which the mean density of the system is 200/9c)- The 
profile at any earlier epoch is then simply a self-similar scaling of this; p/pc being a constant 
function of r/Rmr- 

There is evidence from N-body simulations that the self-similar model breaks down at small 
radii, where some flattening of the density distribution is observed. Following the results of NFW96, 
we therefore flatten the canonical density profile to a logarithmic slope of -1 within a radius 




i/trad, if (2 - dlnA/dluT) > and tcooi < ^mdt 
0, otherwise 



(6) 



2.3 The Cluster Potential 



The gravitational acceleration {g in Equation §) is given by, 

g = -% [Mdarkir) + MicM{r) + Mdep{r)] 



(7) 
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Figure 1: Evolution of the dark matter density in a model with M„jj. = 1.6 x 10^^ ^q- The profile 
evolves in a self-similar way through redshifts 4 (dot-dash), 2 (dotted), 1 (dashed) and (solid). 



of ripRyir, where 7]p is a selectable model parameter. The gas is allowed to evolve within this 
developing dark matter distribution, and the total gravitational potential is computed from the 
combined contributions of dark matter, free gas, and gas deposited by cooling. 

The evolution of the dark matter distribution in physical units, for the case of a cluster with 
virial mass Mmr = 1-6 x 10"*^^ Mq at z = 0, is shown in Figure |l[ In scaled units, p/pc versus 
curves would all fall on top of one another. 

Our approach contrasts with 1-D methods which are based on numerical integration of the 
equations of motion governing the evolution of collisionless dark matter shells (Thoul & Weinberg 
1995). While such methods are more self-consistent than our approach (in the sense that the effect 
of the gas distribution can be included in the calculation of the dark matter dynamics) there is 
a high price to pay in terms of the execution speed of the code, since a large number of dark 
matter shells is required to avoid spurious shocks arising in the gas due to fluctuations in the mass 
distribution. Our approach has the advantage of avoiding this computational overhead whilst 
yielding a good approximation to the expected evolution ~ at large radii the evolution calculated 
by numerical integration will be similar to that given by the analytic solution since the gas and 
dark matter are distributed similarly. Only at relatively small radii, where the gas and dark matter 
distributions can be significantly different, are the two methods likely to differ, but it is precisely 
at this point that 3-D simulations suggest that 1-D treatments break down, presumably due to 
the effect of merging, and the dark matter density profile flattens. We are able to incorporate this 
effect into our model heuristically, by imposing a core on our dark matter distribution. 

Equations [l|-|3| are solved using a standard, first-order, explicit finite-difference Lagrangian 
scheme (Richtmyer &: Morton 1967), with additional terms to account for gravity, cooling, and 
mass deposition. A useful description of the use of a Lagrangian code in relation to cluster cooling 
flows is given by Thomas (1988). Full details of our method are given in Knight (1996). 
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Figure 2: Gas density and velocity profiles from EVOL (solid lines) compared with those predicted 
by self-similar theory (dashed lines; Bertschinger 1985). R^r is the virial radius and Umr = Rvir/t- 
The gas (baryon) fraction was taken to be 0.1. 

3 TESTS OF THE CODE 

The code was subjected to a variety of tests, including comparison with analytically soluble prob- 
lems, previous hydrodynamical simulations, and with the results of another 1-D Lagrangian code, 
kindly supplied by Peter Thomas. Full details of these checks are given in Knight (1996). Here we 
show the results of two of the most important tests - comparisons with the analytical evolution 
derived for a simple 1-D problem by Bertschinger, and with the results of 3-D simulations which 
incorporate the effects of hierarchical merging. 

3.1 Bertschinger Infall Model 

One of the most powerful tests for a hydrodynamical code is the simulation of a self-similar system 
whose behaviour is known analytically. Such a comparison is useful since any departures of the 
simulation from the self-similar behaviour will be attributable to a combination of numerical effects 
and possibly to differences between the model and the simulated system (e.g., different boundary 
conditions). 

Since the evolution of the dark matter profile is given by the infall solution of FG84, in the 
absence of cooling and injection terms the evolution of the gas in the case of a point perturbation 
(n = 0) will follow that derived by Bertschinger (1985). Figure ^ shows the scaled gas density 
and velocity profiles for an EVOL model with n = compared with that predicted analytically. 
Clearly, there is excellent agreement between the profiles. The only significant differences occur in 
the shock region, which is spread out over several shells in the EVOL solutions. 

3.2 Comparison with a 3-D Code 

We have also performed some tests to investigate the differences between EVOL calculations and 
those of a 3-D hydro code (Katz & White 1993). For these simulations of an = 1 universe with a 
CDM initial power spectrum (modelled as an n = — 1 self-similar hierarchy in the case of EVOL), 
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Figure 3: Gas density and temperature profiles calculated by EVOL compared with those calcu- 
lated by a 3-D code (dotted line; Katz & White 1993). 



9% of the mass was taken to be in the form of gas, the rest in dark matter. Radiative cooling 
was included. In the 3-D calculation, a cluster was found to form at the intersection of filaments, 
forming an object with virial radius R^r = 880 kpc and mass Myir = 1.83 x 10^^ Mq by z = 0.13 
when the simulation was halted. The gas density and temperature of this structure is compared 
to that calculated by EVOL in Figure ^ - the initial conditions were taken to be as close to those 
of the 3-D code as possible, and the normalisation of the initial perturbation was determined from 
the virial mass of the final object in the 3-D simulation. 

There is excellent agreement in the shape of the gas density profile at essentially all radii, apart 
from the sudden increase in the density at the centre in the 3-D case due to the central 'cooling 
catastrophe' (which is avoided in EVOL by the distributed gas deposition). The temperature 
profiles agree quite well for r ~ 1-3 Mpc, but that from the 3-D simulation is much more centrally 
peaked than that from EVOL. This is due to the rapid deepening of the potential well as gas near 
the centre clumps under its own self-gravity, becomes even denser and cools even more rapidly, 
leading to the formation of a much deeper cluster potential in the 3-D case. However, Katz and 
White comment that the effect would be reduced if star formation was included (gas can cool and 
stars cannot, reducing the rate at which the potential well can deepen), which is effectively what 
is modelled by the multiphase treatment employed by EVOL. Such central temperature gradients 
are typically absent in 3-D treatments in which cooling is neglected (e.g., Evrard 1990; NFW96), 
most of which yield simulated clusters in which the central regions are essentially isothermal (in 
general agreement with observation, as discussed in Section ^). Given the high sensitivity of the 
results from such treatments to details of the model, we consider the agreement between EVOL 
and 3-D codes to be satisfactory. 

4 THE MODEL GRID 

One of the main advantages of a 1-D treatment is the fast execution time, enabling a much fuller 
exploration of parameter space than would be possible with a 3-D code. This is very useful since 
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the possible parameter space for models of the intraclustcr medium is large, and many of the key 
parameters that can influence the evolution of the ICM are uncertain. In addition, we are especially 
interested in the way that the properties of the ICM scale with system size. We have therefore 
constructed a multidimensional grid that encompasses these parameters, as well as a large range 
in system size, and calculated the evolution of models at a large number (~ 150) of points in the 
parameter space. 

Each model starts with 200 radial shells at 2; = 100, early enough that the assumed initial gas 
profiles have negligible effect on the profiles for z ^ 10. The gas density and velocity were set 
to p(r) = f gas P dark (f)! ^-iid u{r) = HiV (i.e., comoving with the Hubble flow) respectively, where 
Pdarki^) is the initial dark matter density profile (following FG84) and Hi is the Hubble constant 
at the initial epoch. In practice, if a dark matter concentration has formed, then the initial gas 
velocity will actually be somewhat perturbed from the Hubble expansion. However, by starting at 
z = 100, we ensure that this has a negligible effect on the evolution. The initial gas temperature 
was taken to be low (r(r) = 10^ K), ensuring that the initial entropy of the gas is much lower 
than that typically found in the cores of groups and rich clusters. A large initial entropy for the 
gas would effectively introduce a scale size into the evolution, leading to the formation of a gas 
density core. The iron abundance of the gas was set to Z{r) = 0.4, consistent with observed mean 
abundances. In effect, we assume in the present set of models that metals are injected into the 
ICM at some very early stage in cluster evolution - for example, by Population III stars (e.g. Carr 
1994; Ostriker & Gnedin 1996). 

The initial dark matter and gas mass fractions for most of the runs were set to fdark = 0.8 and 
fgas = 0.2 respectively. This represents a baryon fraction of ~ 0.2, higher than that generally 
used for cluster simulations (/& ~ 0.1), and also conflicts (if = 1) with the current restrictive 
constraints on the universal baryon density parameter from primordial nucleosynthesis calculations, 
Qft = 0.05 ± 0.01 /i5Q^ (Walker et al. 1991). However, these constraints are also in severe conflict 
with the observed baryon fractions in clusters of galaxies, 10-20% at 1 Mpc (White &: Fabian 1995), 
rising to ~30% at the virial radius (David, Jones k. Forman 1995). Various possible explanations 
for this difference have been proposed (White et al. 1993), some of which appeal to — 0.3, while 
others suggest that, within an $7 = 1 framework, baryons were preferentially concentrated within 
clusters at an early epoch, as a result of an unknown, non-gravitational mechanism. We take the 
view that, until this issue is resolved, it is preferable to set the baryon fraction in the protoclusters 
to ^ 0.2, rather than produce models in which present day clusters are considerably underdense 
in baryons compared to observation. However, some low baryon density models with fdark = 0.9 
and fgas = 0.1 have also been investigated in order to bracket the observed baryon fraction range 
in clusters and groups. In addition, we have constructed some models to explore the eflfects of a 
lower density Universe (J7o = 0.3) on the evolution. 

By varying the mass enclosed within the virial radius (M^ir) at z = between 7.8 x 10"^^ 
Mq and 1.6 x 10^^ Mq, the grid covers a range of systems from groups with Tyir < 1 keV to 
rich clusters. N-body simulations suggest that the radial infall calculation of the dark matter 
evolution is unphysical at small radii. NFW96 find that the density profile flattens to a logarithmic 
slope of approximately -1 near the centre. Our mass profiles were therefore smoothly flattened to 
p{r) (X within a "core radius" of O.lRyij., in accordance with their result. The results presented 
here are relatively insensitive to the value of r/p, since cooled gas generally becomes the dominant 
gravitational component at small radii. 

The mass deposition parameter was set to ^ = 1.7 for shells where tcool < ctmdi, and ^ = 
otherwise. Within the context of the mass deposition model, values of ^ ~ 1 have both some 
theoretical and observational justification: if ^ is too large the gradient of the specific entropy 
changes sign and the region becomes convectively unstable, effectively reducing ^ back to ~1 
(Fabian, Nulsen & Canizares 1991). The value ^ = 1.7 was chosen since it was found to lead to 
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Parameter Name Units 



Values 



Density parameter Oq — 0.3, 1 

Gas mass fraction fgas — 0.1, 0.2 

Final virial mass M^ir Mq 7.8 x lO^^^ 2.0 x lO^^^ 3 9 x iqI^, 7.8 x lO^^^ 

2.0 X 10^^, 3.9 X 10^"^, 7.8 x 10^"^, 1.6 x 10^^ 

Mass deposition rate ^ — 1, 1.7, 2.5 

Mass deposition threshold Omd — 0.2, 1, 00 

Table 1: Summary of the range of model parameters. 



the formation of simulated deposition profiles of the form M(r) oc r, similar to the form inferred 
from deprojection analyses of cluster cooling flows (Fabian, Nulsen & Canizares 1991). However, 
some simulations with = 1 and ^ = 2.5 have also been performed. 

The mass deposition cutoff at tcooi > otmdt is somewhat arbitrary, although if deposition were 
to occur with ^ ~ 1 out to large radii, the resulting cluster spectra conflict with observations 
- for example, Ginga observations of the Perseus cluster appear to rule out all but a modestly 
multiphase ICM (Allen et al. 1992). 

o/-md — 1 was used for most of the runs, although the effect 
of confining mass deposition to the region tcooi < 0.2t has also been explored. Mass deposition 
reduces the mass of shells within the cooling flow region, helping to improve the resolution of the 
runs ~ shells are dropped when their mass falls below 10^ Mq, and any remaining mass is added 
to a static Eulerian grid which is included in the calculation of the cluster potential. 

5 MODEL EVOLUTION 

The parameter space for the cluster models has many dimensions. It is therefore convenient to 
adopt standard values for the model parameters and then test the effect of alternative values. 
These standard values are taken to be (1) = 1, (2) gas fraction fgas = 0.2, (3) mass deposition 
rate ^ = 1.7, and (4) mass deposition threshold amd = 1. In the following discussion attention is 
focussed on parameters that differ from these standard values. Parameters not explicitly mentioned 
are assumed to have the 'standard' value. 

The evolution of a standard model of mass M^r = 1-6 x 10^^ Mq in physical and "scaled" 
units is shown in Figures § and |5| respectively. The scaled plots were obtained by normalising to 
the characteristic temperature, velocity and entropy: T^r = tJ'mpGMyir/{2kRyir), Umr = Rvir/t, 
and Syir = Tmr/p2^ , where Rmr is the radius within which the mean density, pmr, is 200 times 
the background density at time t. The characteristic temperature, Tyij-, is 8.0 X 10^ K, 3.2 X 10"^ K, 
1.1 X lO'^ K, and 4.4 x 10^ K for models of mass M^ir = 1.6 x 10^^ Mq, 3.9 x 10^^ Mq, 7.8 x 10^^ 
Mq, and 2.0 x 10^^ Mq respectively. Tyir is independent of redshift for n — —1. 

The evolution at large radius is similar to that of the self-similar solution (e.g., Bertschinger 
1985). A particle within the overdense region eventually separates from the Hubble expansion, 
reaching turnaround, u = 0, before collapsing back towards the centre. Before reaching r = 0, 
the fluid particle passes through a shock where the velocity changes discontinuously, after which it 
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Figure 4: Evolution of the gas density, temperature, velocity, entropy, and enclosed gas mass 
fraction of a standard model of mass Myir = 1.6 X 10^^ Mq. Redshifts shown are 0, 0.5, 2, 5 and 
10 (solid, dashed, dotted, dash-dot-dash, and dash-dot-dot-dot-dash lines respectively). 
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Figure 5: Evolution of the scaled gas density, temperature, velocity, entropy, and enclosed gas 
mass fraction for a standard model of mass Myir = 1.6 x 10^^ Mq. Redshifts shown are 0, 0.5, 2, 5 
and 10 (solid, dashed, dotted, dash-dot-dash, and dash-dot-dot-dot-dash lines respectively). The 
thick solid lines in the top panel show the "unsmoothed" dark matter density profile (with cusps) , 
and the profile represented by a power law with an r"^ core applied (see section 2.3). 
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settles to a constant fraction of its radius at turnaround. In contrast, three dimensional simulations 
suggest that the kinetic energy of fluid particles is not completely thermalised after passing through 
the shock. However, the kinetic pressure is only ~10% of the thermal pressure (Evrard 1990) and 
does not affect the results presented here greatly. The position of the shock is marked by a sharp 
increase in gas entropy. Within the shock radius, the entropy rises continuously with radius, 

5 oc r^'^, and convection is suppressed. 

The density profile, while following the dark matter profile quite well at large radius is strongly 
affected by cooling and mass deposition nearer the centre. Gas dropped near the cluster centre is 
replaced by gas of higher entropy from further out, lowering the central gas density. This, coupled 
with the build-up of a large mass of cooled material produces a strong gradient in the enclosed 
gas mass fraction. At z = 0, the gas mass fraction within 0.5 Mpc, 1 Mpc and 3 Mpc is 0.12, 
0.16 and 0.20 respectively. The logarithmic slope of the density profile is found to flatten, from 
approximately -2 to -1, at a radius of ~100 kpc, somewhat less than typically observed core radii 
in rich clusters, Tc — 250 kpc - although many of the systems with the largest core radii may be 
recent merger remnants. As expected, the effects of cooling are most dramatic at early epochs, 
when the gas density is higher and the effect of mass deposition consequently more widespread, 
leading to a reduction in the gas mass fraction throughout the cluster. At later times, cooling is 
only important near the cluster centre. 

Figure ^ shows that the temperature generally declines continuously with increasing radius, 
even within the cooling flow region. This is not unexpected since the "hot phase" of such a flow 
is effectively being modelled (the cooling gas is dropped from the ICM at each timestep), and the 
gas will be adiabatically compressed as it flows inward. Recent studies of multiphase cooling flows 
within mass profiles with small cores confirm that the temperature of the hot phase will remain 
approximately constant within the cooling region (Waxman &; Miralda-Escude 1995), while if the 
gravitating matter distribution is very sharply peaked, as found here, the temperature is expected 
to rise (White & Sarazin 1987). Beyond this region, the temperature gradient steepens with radius 
until the sudden drop at the shock. Similar temperature profiles have been found in a number of 
three dimensional hydrodynamical simulations (e.g. NFW96). 

Observations of the temperature structure of clusters have been mostly restricted to the nearest 
and brightest systems until recently. Observations with coded mask and collimated instruments 
indicate that the Coma cluster (Hughes, Gorenstein & Fabricant 1988; Watt et al. 1992) has 
an approximately isothermal core within '^l Mpc of the cluster centre with a steep temperature 
decline at large radius, whilst the Perseus cluster (Eyles et al. 1991) shows a modest temperature 
decline outside the cooling region, with the emission well modelled by a power-law temperature 
profile, T(r) oc r-0-30±o.08_ j^ecent results from ROSAT and ASCA (e.g. Allen, Fabian & Kneib 
1996; Markevitch & Vikhlinin 1996; Loewenstein & Mushotzky 1996) suggest that most reasonably 
relaxed clusters have fairly flat temperature distributions within the central ~1 Mpc (apart from 
the common presence of a central cooling flow), although steep gradients at larger radii appear to 
be present in at least some high temperature clusters (Markevitch 1996). 

Figure |6| shows the final {z = 0) scaled gas profiles for different mass models. The potentials 
of these clusters have evolved in a self-similar fashion, as described in section 2.3. If the gas has 
followed suit, then the profiles in Figure |6| should coincide. At large radius, where there is no 
additional physics to break the scaling, the scaled profiles are very similar to each other. However, 
significant differences occur near the cluster centre, where cooling is important. 

6 THE COOLING FLOW 
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Figure 6: Final scaled gas density, temperature, velocity, entropy, and enclosed gas mass fraction 
profiles for standard models of mass M^ir = 2.0 x lO^^, 7.8 x 10^^ 3.9 x 10^^ and 1.6 x 10^^ Mq 
(dot-dash, dotted, dashed, and solid line styles respectively). 
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6.1 Mass Deposition Profile 

As discussed in Section the mass deposition parameter ^ was chosen so that the mass deposition 
profile is similar to observations, M(r) oc r (Fabian, Nulsen & Canizares 1991). The final cumula- 
tive mass deposition rate is shown plotted against radius for the Myir = 1.6 x 10"*^^ Mq model in 
Figure 0a - within the mass deposition region M(r) oc r, as expected. The total mass deposition 
rate for this model is higher than generally found observationally - M ~ 2000 Mq compared to 
typical observed values of < 1000 Mq yr~^, although there are a few exceptional cases with much 
higher inferred mass deposition rates (e.g. Schindler et al. 1996). The density profile of the cooled 
material is centrally peaked and constitutes the dominant gravitational component at small radii 
(see Figure HI). The cluster potential is therefore sufficiently peaked that the cooling flow becomes 
gravity dominated near the cluster centre, and the temperature of the hot phase increases due to 
adiabatic compression as it flows inwards. As a result of this dominant cooled gas component in 
the core, the circular velocities near the cluster centre are Vdrc — 1500 km s~^ in high mass models 
- much larger than the velocity dispersions ( < 400 km s~^, Oegerle &: Hoessel 1991) observed in 
the envelopes of the cD galaxies which are commonly found at the focus of cluster cooling flows. 
This discrepancy could be reduced by mechanisms not incorporated in the present models, such 
as cluster mergers distributing the material throughout the core, or preheating of the cluster gas, 
which reduces the mass of cooled material (see Paper II). 

Several possible mechanisms for lowering the mass deposition rate have been explored: lowering 
the baryon fraction, fgas, artificially reducing the radiative cooling rate, and varying the model 
parameters governing the deposition of gas in the cooling fiow (^ and Omd)- The final {z = 0) mass 
deposition profiles for some of these models are shown in Figure |^-d. Except for the models with 
low baryon fraction, varying the above parameters has little effect on the total mass deposition 
rate. To some extent, there is a feedback mechanism which keeps the mass deposition rate at an 
approximately constant value, regardless of the values specified for the cooling flow parameters - if 
the cooling rate is reduced for example, the gas simply flows further into the cluster centre (where 
the density is higher) before dropping out of the flow. The models with reduced baryon fraction 
have significantly lower mass deposition rates, but the gas mass fraction in these cases is in confiict 
with observations (White & Fabian 1995), even if the ICM is multiphase (Gunn & Thomas 1995). 

Since our attempts at reducing the mass deposition rate by varying the parameters of the 
code have been unsuccessful, we conclude that the most likely explanations for the discrepancy 
between the mass deposition rate of the models and values inferred from observations are related 
to limitations of the model. For example, (1) the spherical geometry assumption, and the neglect 
of angular momentum that would inhibit gas flowing into the cluster centre, (2) the lack of cluster 
mergers that would periodically disrupt the cooling fiow, and (3) turbulence and thermal conduc- 
tion within the core. Since both angular momentum and merging would tend to reduce the central 
density, and hence the mass deposition rate, we believe that our models represent 'maximal' cooling 
fiows, giving a useful upper bound to mass deposition rates that might be found observationally. 
This view receives some support from the observed variation of mass deposition rate with cluster 
temperature, as discussed in the next section. 

The high mass deposition rates of the standard models have a knock-on effect on other quan- 
tities, such as the emission-weighted temperature, and X-ray luminosity. Since the lower mass 
deposition rates of real clusters are likely to be due to physical processes not modelled here, we 
have shown the effect, in some of the succeeding plots, of arbitrarily reducing these mass deposition 
rates by a factor of 5. 
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Figure 7: Final cumulative mass deposition rate for the My^r = 1-6 x 10^^ Mq model. Panels 
show (a) ^ = 1.7, amd = 1-0, fgas = 0.2, (b) ^ = 1.7, 1.0 and 2.5 (solid, dashed and dotted lines), 
(c) ctmd = 1-0) 0-2 and 5.0 (solid, dashed and dotted lines), and (d) fgas = 0.2 and 0.1 (solid and 
dashed lines respectively). 
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Figure 8: Final cooled gas, dark matter and gas density profiles for the M„, 
model (solid, dashed and dotted lines respectively). 
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6.2 The Effects of Merging 

Figure ^ shows the variation of total mass deposition rate with emission- weighted temperature for 
models of different mass. Two effects are apparent: (i) the mass deposition rates of the models are 
higher than those found observationally, and (ii) although the deposition rate increases with mass 
in the models, the trend is less steep than in observed clusters. For high mass systems, the models 
shown in Figure ^ appear to represent an upper bound to the observations, suggesting that some 
additional variable (perhaps cluster mergers) interrupts the idealised behaviour of the simulations. 
For example, if cooling flows are disrupted by major mergers and the time between such mergers 
is comparable to the timescale for the cooling flow to re-establish itself following a merger, then a 
wide range of mass deposition rates for clusters of a given richness might be expected. 

In order to test this possibility, and also to explore the impact on our other results, we have 
performed a series of simulations in which the cooling flow is disrupted by raising the entropy 
of the gas within the cluster core (r < O.lR^r)- The effect of this is to force gas out of the 
cluster core, temporarily reducing the mass deposition rate - Figure IC shows the evolution of 
the mass deposition rate for a simulated rich cluster in which the cooling flow is disrupted several 
times during its evolution using this procedure. To determine whether cooling flow disruption 
can account for the scatter in the variation of mass deposition rate with cluster temperature, we 
have constructed some models in which the cooling flows are disrupted periodically, as would be 
expected in a Universe in which structures grow largely by hierarchical merging. We assume that 
clusters typically experience a major merger each time the Universe doubles its age (see e.g. Lacey 
&: Cole 1993; Kauffmann &; White 1993), but we also introduce some scatter in the merger times, 
such that the nth major merger occurs on average at time t", where = 2t^~^, with the actual 
merger time being chosen randomly from the range {t"'~^ + t^)/2 (t" -|- t^~^^)/2. The results 
from 80 such models are shown in Figure 0. 
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Figure 9: Final total mass deposition rate plotted against emission-weighted temperature for 
models of mass Myir = 7.8 x 10^^ 2.O x lO^^ j^^^ 3 9 ^ iqIS j^^^ 7 3 ^ iqIS M©, 2.0 x 10^^ 
M0, 3.9 X 10^^ Mq, 7.8 X 10^^ Mq and 1.6 x 10^^ Mq. The lower line shows the effect of scaling 
the mass deposition rate down by a factor of five. The observational data are taken from White, 
Jones & Forman (1996). 
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Time (Gyr) 

Figure 10: Evolution of the total mass deposition rate for a model of mass Myir = 1.6 XlO^^ Mq m 
which the cooling flow is disrupted using the procedure described in the text at several redshifts. 
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Figure 11: Final total mass deposition rate plotted against emission- weighted temperature for the 
cooling flow disruption models described in the text. The filled circles are observational data taken 
from White, Jones &z Forman (1996), and the arrows indicate model clusters with M < 1 Mq yr~^. 
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Mean cluster 


Final virial mass 


Final mean model 




temperature^ 


of model 


temperature 




T/keV 




T/keV 


Hickson 62 


0.96 


7.8 X 10^^ 


1.1 


A1060 


2.55 


2.0 X lO^-^ 


2.2 


MKW 3 


3.0 


3.9 X lO^'' 


3.4 


A1795 


5.34 


7.8 X 10^^ 


5.4 


1 Ponman & 


: Bertram 1993; 


Yamashita 1992 





Table 2: Mean temperatures for the clusters in Figure 12a-d. 



Clearly, disruption of the cooling flow can substantially reduce the present day mass deposition 
rates for the more massive systems since, for these systems, the timescale for the cooling flow to 
re-establish itself is comparable to the time between mergers. However, the low mass systems have 
substantially lower cooling times in their cores and the cooling flow is able to re-establish itself 
much faster, hence they are, on average, little affected by cooling flow disruption. For massive 
systems, the scatter in the mass deposition rate for a given temperature appears to be similar to 
that observed, but high resolution 3-D simulations with radiative cooling are needed to definitively 
establish the effects of cluster merging on cooling flows. 



6.3 Evolution of the Cooling Flow 



The evolution of the mass deposition rate is shown for a range of system masses in Figure 12. The 
mass deposition rate is substantially (a factor of 2-3) larger at z = 5 than at z = 0, decreasing 
with time due to the inflow of higher entropy gas from larger radii to replace that deposited at the 
cluster centre. This lowers the central gas density, which in turn reduces the mass deposition rate. 
As a result, most deposition occurs at high redshift, with ~ 55% of the total amount dropped by 
z ~ 1. 



7 THE GAS DISTRIBUTION 

7.1 The X-ray Surface Brightness Profile 



Final z = bolometric surface brightness profiles for different mass models are shown in Figure |1^, 
together with observed profiles for groups and clusters of comparable mass (Table |2|) . The obser- 
vational profiles are based on ROSAT PSPC data converted to bolometric flux using ROSAT or 
Gm^a temperatures, assuming isothermal gas (Cannon, Ponman & Navarro, in preparation). The 
model profiles were calculated by integrating the emission from the shells along each line of sight, 
using the cooling function of Raymond, Cox & Smith (1976). An additional component due to 
gas cooling to zero temperature at a rate proportional to the local mass deposition rate was also 
included. 

The surface brightness profiles of the more massive systems are similar to, but slightly steeper 
than, those of observed clusters of similar mass: at the centre, the very large contribution of mass 
deposition to the cluster luminosity is clearly visible as a region of excess emission. Beyond this 
region, the surface brightness profile gradually steepens with radius. The shape of the EVOL 
surface brightness profiles of the richer systems are quite similar to the observed profiles outside 
the central region. However, unlike the EVOL profiles, there is no evidence for a region of greatly 
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Figure 12: Evohition of the mass deposition rate for models with (from top to bottom) My, 
1.6 X 10^^ Mq, 3.9 X lO^'' Mq, 7.8 x 10^^ Mq, and 2.0 x 10^^ 
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Figure 13: Final bolometric surface brightness profiles for systems of mass (a) Myj,. = 7.8 x 10 
M©, (b) M^ir = 3.9 X 10^^ Mq, (c) Myir = 2.0 x 10^^ M© and (d) 7.8 x lO^^ M©. The dotted 
line shows the surface brightness when the mass deposition luminosity is ignored. The circles 
indicate the surface brightness profiles of the following clusters observed by the ROSATVSPG (D. 
Cannon, private communication; Ponman & Bertram 1993): (a) A1795, (b) MKW3, (c) A1060, 
and (d) Hickson 62. Note that the models have not been optimised to match the data, which are 
illustrative only. The relationship S oc is shown for reference, and Hq = 50 km s~^ Mpc""^ is 
assumed. 
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enhanced emission due to mass deposition within the coohng flow near the centre of observed clus- 
ters (though the isothermal assumption may in practice lead to an underestimate of the bolometric 
surface brightness by ~ 20% within the cooling flow region). The mismatch between the models 
and data near the centre persists even if the numerical parameters in EVOL are adjusted to conflne 
the mass deposition to a region closer to the cluster centre. The discrepancy could be resolved if 
the cooling flows of real clusters are periodically disrupted (Section ^). As can be seen from the 
flgure, the model proflles of the massive systems match the data much better if the emission from 
the model cooling flow is omitted. 
A fit of the function 

-3/3^,t+l/2 

(8) 

to the underlying cluster emission (i.e., excluding the additional luminosity due to mass deposition) 
between r = and O.bRvir yields a Pfu value of ~0.72. This is slightly larger than inferred from 
observations of rich clusters, which generally yield <(3fit>~ 0.67 and hence Sx oc r~'^ outside the 
core (Jones & Forman 1984). 

For lower mass clusters, the emission from the models falls above that found observationally, 
with the largest discrepancy being found in the smallest systems. This suggests that the gas 
fraction in the inner regions declines from rich clusters to groups, as found by David, Jones & 
Forman (1995). 
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7.2 (3fit vs Temperature and the Gas Entropy 

Figure |l^ shows the variation of Pfu with emission- weighted temperature at z = for models of 
different mass. Also shown are values of Pfu from an observed sample (David, Jones &: Forman 
1995). 

The results confirm those of Section 7.1 - there is a modest trend for a decrease in f3fit with 
decreasing T, although /3/it is found to drop substantially for the lowest mass system, probably 
due to the importance of line cooling. Reasonable agreement with observations is found for the 
rich systems which have iSfu — 0.68, very similar to the values found observationally (/3/jt — 0.67). 
There is, however, a large discrepency between the simulated clusters and observations for the 
poorer clusters. Observationally, Pju is found to decline steeply with decreasing T, which suggests 
the presence of an additional mechanism operating in poorer systems that has not been included 
in the simulations presented here. 

Another pointer to the need for additional physical processes in low mass systems is the value of 
the entropy of the intracluster gas. Gas in the inner regions of X-ray bright galaxy groups (outside 
any cooling fiow) is observed to have an entropy corresponding to S = T/ns ~ 200 — 300 keV cm^ 
(Bourner and Ponman, paper in preparation) - very similar to that found in cluster cores. In 
the EVOL models presented here, the gas entropy at the corresponding point in rich systems has 
S 200—400 keV cm^, as observed, but in the low mass models (T 1 keV) it has S* < 50 keV cm^. 



7.3 The Gas Mass Fraction Profile 



Figure 15 shows the final enclosed gas mass fraction proflles for different mass models. In each case, 
the processes of cooling and mass deposition, particularly at early epochs, lead to the development 
of a profile that rises continually with radius (gas which has cooled out is counted as mass, but 
not as 'gas'). Such a rise is similar to that inferred by White & Fabian (1995) from a sample of 
Abell clusters, and also agrees reasonably well with the results of David, Jones & Forman (1995). 
The mechanism behind this is interesting - since the entropy of the gas increases with radius, mass 
deposition near the cluster centre leads to the inflow of higher entropy gas from further out in 
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Figure 14: Final f3fit plotted against emission- weighted temperature for models of mass Myir = 
2.0 X 10^''^ Mq, 3.9 X 10^3 Mq, 7.8 x lO^^ M©, 2.0 x 10^^ Mq, 3.9 x 10^^ Mq, 7.8 x 10^^ M© and 
1.6 x 10^^ Mq. The observational data (crosses) are taken from David, Jones & Forman (1995). 



25 








J I L_ 



J I L 



J I L 



-2 -1 

Log(r/R^.^) 







Figure 15: Final enclosed gas mass fraction profiles for models of mass Myir = 1.6 X 10^^ Mq, 
3.9 X 10^^ Mq, 7.8 X 10^^ M© and 2.0 x 10^^ Mq (solid, dashed, dotted, dot-dash lines respectively). 
The horizontal dotted line shows the global value. 
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the cluster which, in turn, lowers the central gas density producing a gas "core" . The total mass 
density profile is, however, quite sharply peaked despite the core in the dark matter profile, being 
dominated by "cooled" gas at small radii (Figure ^). The net result is a gas mass fraction profile 
that increases with radius. The slight decline in fgas seen outside Rmr arises from the fact that 
the gas density drops sharply at an outward propagating shock, whilst the dark matter density 
(which is represented by a cored power law model, as discussed in section 2.3) declines smoothly 
to the background value. 

Although central cooling and consequent gas inflow does produce a rising gas fraction profile, 
in accord with observations, it involves mass deposition rates, particularly at large redshift, which 
are much larger than those typically inferred for nearby clusters. Alternatively, there are several 
possible mechanisms, such as major mergers, and energy injection from the cluster galaxies, which 
may be effective at forcing gas out of the cluster core - this would reduce the gas fraction in the 
inner regions, and should also lower the mass deposition rate in the cooling flow. 

Figure 15 also shows that fgas in the inner regions declines somewhat with decreasing system 
size, due to the greater importance of cooling in lower mass systems. 



7.4 The Effect of a Fully Multiphase ICM 

Since the cooling time threshold for the onset of multiphase mass deposition is somewhat arbitrary, 
we have examined the alternative possibility that the bulk of the ICM is multiphase by constructing 
some models in which mass deposition is allowed to occur throughout the cluster {omd = oo). The 
final (z = 0) surface brightness profile for such a model of mass M^ir = 7.8 X 10^^ Mq is shown 
in Figure The shape of the surface brightness profile is similar to that of observed clusters, 
but the contribution of the cooling gas to the cluster luminosity is appreciable at all radii (~ 70% 
of the total luminosity). This reduces the gas fraction required for a given cluster luminosity by 
~ 70% and would be a possible mechanism for alleviating the "baryon catastrophe" associated 
with clusters of galaxies (Gunn & Thomas 1995). However, it was shown some time ago by Allen 
et al. (1992) that the spectral properties of the intracluster emission in the Perseus cluster are 
inconsistent with such large scale multiphase structure, and this conclusion has been confirmed for 
other systems by ASCA (Mushotzky et al. 1996). 



8 X-RAY LUMINOSITY 

8.1 Luminosity- Temperature Relation 



Figure [17| shows the bolometric luminosity, Lf,oi, plotted against emission- weighted temperature, 
Tijoi, for various models. The calculation of T^o; includes the contribution of mass deposition within 
the cooling flow. For the hottest models, this "cooling luminosity" lowers the emission- weighted 
temperature by up to ~30%. If the structure of all clusters is self-similar then scaling relations 
can be derived for the evolution of mean cluster properties. One of the predictions of self-similarity 
is a relationship between bolometric luminosity and temperature of the form L^o; cx (NFW96), 
much flatter than the observed relationship, oc 7"2.8i±o.i8 ^white, Jones Sz Forman 1996). 

The slope of the L : T relation for the infall models considered here is 2.1 (2.2 if the mass 
deposition rate is reduced by a factor of five) , which is still flatter than that found observationally. 
In order to investigate the effects of cooling on the lower mass systems, some simulations were 
performed in which radiative cooling was turned off during the cluster evolution. The results of 
viewing the X-ray emission from such a model at z = are shown in Figure |l8|. (The surface 
brightness profiles of these models are so steep that the luminosity diverges at small radii so, for 
these models only, emission within a radius of 0.05Rmr has been ignored.) The greater importance 
of line cooling to the bolometric luminosity for the low mass systems causes the L : T relation to 
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Figure 16: Final bolomctric surface brightness profiles for a fully multiphase ICM model of mass 
7.8 X 10^'^ Mq. The dotted line shows the surface brightness when the mass deposition luminosity 
is ignored. The circles indicate the surface brightness profile of A1795 observed by the ROSAT 
PSPC (D. Cannon, private communication). 
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Figure 17: Bolomctric luminosity plotted against temperature for models of mass Myir = 7.8 x 10^^ 
Mq, 2.0 X 10^3 Mq, 3.9 X 10^"^ Mq, 7.8 x 10^^ Mq, 2.0 x 10^^ Mq, 3.9 x 10^^ M©, 7.8 x 10^^ 
Mq and 1.6 X 10^^ Mq. The lower line shows the effect of scaling the mass deposition rate down 
by a factor of 5. The data points are from Ponman et al. (1996), and the dashed line (which has 
arbitrary normalisation) shows the relation, L (xT'^, expected for self-similar cluster structure. 
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Figure 18: Bolometric luminosity plotted against temperature for models in which radiative cooling 
has been turned off during the evolution. The masses of the clusters arc 

2.0 X 10^3 Mq, 3.9 X 10^^ Mq, 7.8 x lO^^ Mq, 2.0 x 10^^ Mq, 3.9 x 10^"^ Mq, 7.8 x 10^^ Mq 
and 1.6 X 10^^ Mq. The lower line shows the effect of setting the metal abundance to zero when 
calculating the luminosity and mean temperature. The data points are from Ponman et al. (1996), 
and the dashed line (which has arbitrary normalisation) is the prediction for self-similar cluster 
structure, L oaT"^. 
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flatten at low temperatures. However, when line radiation is neglected, the logarithmic slope of the 
L : T relation is 1.95, very close to the value of 2 expected if the cluster population is self-similar, 
with different mass clusters having profiles of the same form. 

This test shows that the fact that the slope of the L : T relation in Figure 17 is so similar to 
that expected for self-similar scaling is somewhat fortuitous. It occurs because the flattening in 
the relation due to the increased contribution of line emission and mass deposition in low mass 
systems is compensated by the fact that this extra mass deposition lowers the gas fraction, which 
steepens the L : T relation. The net result is an L : T index only slightly steeper than 2. 

The second effect which is apparent in Figure 17 is that the normalisation of the luminosity of 
the simulated clusters is larger than that observed. This discrepancy is partly due to the fact that 
the mass deposition rates are considerably larger than typically inferred from observations. The 
contribution of this cooling material to the overall cluster luminosity is signiflcant (see Figure 13) 
- this is a generic feature of the simple infall models, possibly compounded by the simplified radial 
cooling flow. 

In Figure |l^ we show the L : T relation for the set of merger models discussed in Section |6.2| . 
The effects of merging are to substantially reduce the luminosity of the more massive (T > 3 keV) 
systems, bringing them into line with observed clusters. However, as noted in Section |6.2| , the 
effect of merging on lower mass systems is much less marked, due to the speed with which the 
system recovers from a merger, and re-establishes its cooling flow. 



8.2 Luminosity Evolution 



The evolution of the bolometric luminosity for models with a range of masses is shown in Figure 2C . 
X-ray emission is generally detectable out to 0.2-0.5 Mpc and ~2 Mpc for groups and clusters 
respectively (Mulchaey et al. 1996; Ponman et al. 1996; White &: Fabian 1995), corresponding to 
a radius of ~ 0.5Rmr- Accordingly, the luminosities given here were calculated by summing the 
emission from shells within a radius of 0.5Rmr- A line showing the evolution expected if the gas 
distribution is self-similar (e.g., if Pgasii") oc Pdark{r)) is also shown in Figure 

L,o/oc(l + z)(i3+^")/('^+2"). (9) 

This relation is found to give a very good match to the evolution of the models with no radiative 
cooling. However, as shown in the figure, models with radiative cooling evolve much less strongly 
- the cluster luminosity rises by ~50% out io z = 0.5, while Equation |9| gives a rise of ~80% for 
n = — 1 over the same interval. The difference is predominantly due to evolution of the central 
gas mass fraction: at earlier epochs, the characteristic density of the virialised region is higher and 
mass deposition more widespread, resulting in a lower gas mass fraction within any fixed fraction 
of the virial radius (see Figure ^). As a result, the gas mass fraction rises with time compared 
to the self-similar scaling, reducing the evolution of the luminosity. Such modest (~30% out to 
z = 0.3) positive evolution of the X-ray luminosity of rich clusters is consistent with recent results 
from clusters observed during the ROSAT A\\ Sky Survey (H. Ebeling, private communication). 



9 THE EFFECT OY VLq < I 

In order to test the extent to which the results are sensitive to the assumed value of ^q, a few 
models have been considered in which < 1. In this case, the evolution of the dark matter density 
departs significantly from the self-similar behaviour at redshifts z < Zcrit = ^ — 1 (Lacey &: Cole 
1993). To simulate the approximate evolution of a cluster in an open Universe, it is assumed that 
the dark matter density behaves exactly as in the Vt = 1 case for z > Zcrit-, whilst for z < Zcrit 
the dark matter density beyond a radius of rcru = 5Rmr (roughly the radius at which material is 
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Figure 19: Final bolometric luminosity plotted against emission-weighted temperature for cooling 
flow disruption models. The data points are from Ponman et al. (1996), and the dashed line (which 
has arbitrary normalisation) shows the relation, L ocT^, expected for self-similar cluster structure. 
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Figure 20: Evolution of the bolometric luminosity for models of mass (from top to bottom) M^r = 
1.6 X 10^5 Mq, 3.9 X 10^'' Mq, 7.8 x 10^^ Mq, and 2.0 x 10^^ Mq. The dashed line shows the 
self-similar prediction. 
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just turning around from the Hubble expansion) is constant, and equal to the mean density of the 
universe at that epoch. 

The treatment of the dark matter within rcrit after z = Zcrit is somewhat problematic, due to 
uncertainty in the evolution of the dark matter "core". We consider two cases: (1) at z = Zcrit, 
the profile within rcrit is frozen and there is no subsequent evolution, and (2) the dark matter 
profile within rcrit is allowed to evolve in the same way as in the critical density case. Clearly the 
behaviour of the low models in both cases will be unrealistic at large radius - the density at rcrit 
will suddenly drop into the background, whereas in reality one would expect a gradual steepening 
of the density profile at large radius (Hoffman & Shaham 1985; Crone, Evrard & Richstone 1994). 
However, the dynamics in the central region, which contributes most of the luminosity, will be 
relatively unaffected by this. 

Figure 21 compares the evolution of the luminosity for two ^Iq = 0.3 models (which differ in the 
handling of the core evolution as discussed above) and an = 1 model, all with M^r ~ 1-0 x 10^^ 
Mq at z = 0, and a gas fraction of 0.3. Clearly the dependence of the luminosity evolution on Oq, 
and on the treatment of the core in the Qq = 0.3 models, is slight. The weakness of the dependence 
on r^o comes about because the luminosity is dominated by gas in the central regions, which has 
a very similar profile in the = 0.3 and 0,q = 1 cases. Significant differences occur at large radii, 
but these regions contribute little to the X-ray flux. 



10 CONCLUSIONS 

The 1-D simulations presented above, including an evolving potential and the effects of radiative 
cooling and distributed mass deposition, are reasonably successful at reproducing many of the 
observed properties of galaxy clusters. Temperature and density profiles are broadly realistic, and 
the inflow of higher entropy gas into the cluster core as a result of radiative cooling, leads to a gas 
fraction which rises with radius, in accord with observations. It also results in X-ray luminosity 
evolution which is milder than the self-similar prediction. The strong positive evolution of the 
latter is known to conflict with observation (Castander et al. 1994). 

There are, however, some significant discrepancies with observation. The cooling flows pro- 
duced by EVOL simulations are much stronger than those typically observed. This result is robust 
to changes in the values of numerical parameters in the simulations. The effects of subcluster 
merging, which are believed to periodically disrupt cooling flows, and are of course absent from a 
1-D code, are likely to account for the lower cooling rates (and the substantial scatter for clusters 
of a given X-ray luminosity) seen in massive clusters. However in lower mass systems, the cooling 
flow re-establishes itself very quickly after disruption. 

A number of further disagreements with observation for low mass systems suggest the need 
for additional physics in the models. The simulations reported here do give somewhat flatter 
surface brightness profiles (low fj values - Figure 14) and lower gas fractions (Figure ^) in low 



mass systems, but the effects are much weaker than those observed. Also the slope of the L : T 



relation does not steepen below T ~1 keV in the way which is observed (Figure 17). Perhaps most 
suggestive of all, is the observation that the simple models described in this paper underpredict 
the entropy of the gas in groups and poor clusters. 

These differences suggest that some additional effects, which are most significant in low mass 
systems, have increased the entropy and reduced the density of gas in the inner regions of low mass 
systems. Possible mechanisms include heating of the gas to a temperature > 10^ K, either as the 
result of an early preheating era (Couchman &: Rees 1986; Meiksin & Madau 1993) or the injection 
of galaxy winds (White 1991; Metzler & Evrard 1994), or increased efficiency of star formation in 
groups (David, Jones & Forman 1995). 

The role of these additional mechanisms is explored in Paper H. 



34 



m 

CD 
O 



o 



25 



20 



15 



10 - 



5 








4 







Redshift 



Figure 21: Evolution of the bolometric luminosity for models with (1) J7o = 1 and M^r = 1-6 x 10^^ 
Mq (solid line), (2) J7o = 0.3, M^j^ = 8-6 x 10"^^ -^0, and dark matter profile within frozen for 
z > Zcrit (dashed line), and (3) J^o = 0.3, Alvir = 1-3 x 10^^ Mq, and dark matter profile within 
allowed to evolve as in the = 1 case (dotted line). 
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